# Load Packages -----------------------------------------------------------
library(psych)
library(stargazer)

# Load Data ---------------------------------------------------------------
# Set file path
setwd("")

load("drop_in_ocean_data.rdata")

# PCAs ------------------------------------------------------------
# define variables
direct.vars <- c("pol.stop", "court.appear", "prison.time")
reform.vars <- c("felon_benefits", "felon_voting", "cit_osight", "blm_support", "pub_defend")
eval.vars <- c("p.race.fair", "p.exces.force", "p.account", "court.fair.all")

# estimate pca for IV and store values
dir.pca <- pca(subset(cjs.df, select = direct.vars), missing = T)
cjs.df$direct_exp_pca <- as.numeric(dir.pca$scores)
cjs.df$direct_exp_pca <- (cjs.df$direct_exp_pca + abs(min(cjs.df$direct_exp_pca, na.rm = T)))/(max((cjs.df$direct_exp_pca + abs(min(cjs.df$direct_exp_pca, na.rm = T))), na.rm = T))

# estimate pca for DV and store values
cjs.pca <- pca(subset(cjs.df, select = c(eval.vars, reform.vars)), missing = T)
cjs.df$cjs_att_pca <- as.numeric(cjs.pca$scores)
cjs.df$cjs_att_pca <- (cjs.df$cjs_att_pca + abs(min(cjs.df$cjs_att_pca, na.rm = T)))/(max((cjs.df$cjs_att_pca + abs(min(cjs.df$cjs_att_pca, na.rm = T))), na.rm = T))


# Race Analyses -----------------------------------------------------------
# estimate model
cjs_att_mod_race <- lm(cjs_att_pca ~ direct_exp_pca*black + rr_sc*black + 
                           peer.felony*black + age_sc*black + man*black + 
                           educ_sc*black + inc_sc*black + pid7_dem*black,
                         data = cjs.df, weights = wts_bw_stk)
summary(cjs_att_mod_race)


# By Gender (Black) -----------------------------------------------------
# Estimate Model
cjs_att_mod_gender_blk <- lm(cjs_att_pca ~ direct_exp_pca*man + rr_sc*man +
                           peer.felony*man + age_sc*man + educ_sc*man + inc_sc*man + pid7_dem*man,
                         data = cjs.df, subset = black == 1,  weights = wts_bw_stk)
summary(cjs_att_mod_gender_blk)


# By Gender (White) -----------------------------------------------------
# Estimate model
cjs_att_mod_gender_wht <- lm(cjs_att_pca ~ direct_exp_pca*man + rr_sc*man +
                           peer.felony*man + age_sc*man + educ_sc*man + inc_sc*man + pid7_dem*man,
                         data = cjs.df, subset = black == 0,  weights = wts_bw_stk)
summary(cjs_att_mod_gender_wht)


# Tables -------------------------------------------------------
stargazer(cjs_att_mod_race, cjs_att_mod_gender_blk, cjs_att_mod_gender_wht,
          title = "Direct Experiences Effect Carceral State Attitudes by Race and Gender, Principal Components Model",
          dep.var.caption = "",
          dep.var.labels.include = FALSE,
          column.labels = c("", "Black", "White"),
          covariate.labels = c("Direct Experiences (PCA)", "Experiences (PCA)*Black", "Experiences (PCA)*Man",
                               "Black", "Man", "Black*Man",
                               "Racial Resentment", "Resentment*Black", "Resentment*Man",
                               "Peers with Felony Convictions", "Peers*Black", "Peers*Man",
                               "Age", "Peers*Black", "Peers*Man",
                               "Education", "Education*Black", "Education*Man",
                               "Income", "Income*Black", "Income*Man",
                               "Partisanship", "Partisanship*Black", "Partisanship*Man"),
          order = c(1, 10, 18,  2, 6, 14, 3, 11, 19, 4, 12, 20, 5, 13, 21, 7, 15, 22, 8, 16, 23),
          model.numbers = F,
          no.space = T, header = F,
          notes = c("OLS regression results with standard errors in parentheses. Analyses weighted. Measures scaled 0-1."), 
          notes.align = "l", 
          intercept.bottom = T, initial.zero = F, label = c("pca_out"),
          type = "text",
          # out = c("tableA6.tex"),
          digits = 3, omit.stat = c("f", "adj.rsq"), df = F,
          align = T, star.char = c("*"), star.cutoffs = c(0.05))

